Astronomy & Astrophysics manuscript no. 


February 2, 2008 


(DOI: will be inserted by hand later) 





Temporal evolution of magnetic molecular shocks 

II. Analytics of the steady state and 
semi-analytical construction of intermediate ages. 

P. Lesaffre 1 ' 2,3 , J.-P. Chieze 3 , S. Cabrit 4 , and G. Pineau des Forets 5 ' 6 

■ 1 Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK 

2 University of Oxford, Department of Astrophysics, Oxford 0X1 3RH, UK 
CN ' 3 CEA/DAPNIA/SAp Orme des Merisiers, F-91191 Gif-sur-Yvette Cedex 

4 LERMA, UMR 8112 du CNRS, Observatoire de Paris, 61 Av. de l'Observatoire, F-75014 Paris 

5 IAS, UMR-8617 du CNRS, Universite Paris-Sud, bat. 121, F-91405 Orsay 

6 LUTH, UMR-8102 du CNRS, Observatoire de Paris, F-92190 Meudon Cedex 

oo i 

! Received December 15, 2003; Accepted June 10, 2004 



> 



Abstract. In the first paper of this series (Paper I) we computed time dependent simulations of multifluid shocks 
qq 1 with chemistry and a transverse magnetic field frozen in the ions, using an adaptive moving grid. 

qq ' In this paper, we present new analytical results on steady-state molecular shocks. Relationships between density 

iy~) , and pressure in the neutral fluid are derived for the cold magnetic precursor, hot magnetic precursor, adiabatic 

£ — ' shock front, and the following cooling layer. The compression ratio and temperature behind a fully dissociative 

| adiabatic shock is also derived. 

To prove that these results may even hold for intermediate ages, we design a test to locally characterise the validity 
of the steady state equations in a time-dependent shock simulation. Applying this tool to the results of Paper I, 
we show that most of these shocks (all the stable ones) are indeed in a quasi-steady state at all times, i.e. : a 
given snapshot is composed of one or more truncated steady shock. Finally, we use this property to produce a 
construction method of any intermediate time of low velocity shocks (it < 20 kms -1 ) with only a steady-state 
code. In particular, this method allows one to predict the occurrence of steady CJ-type shocks more accurately 
than previously proposed criteria. 
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C$ ' 1. Introduction terpret observed spectra. Nevertheless, steady-state codes 

■ 1 I 1 have their limits. Observed magnetic molecular shocks 

In a previous paper jLesaffre et al\ |200i hereafter may nQt be fuUy in steady _ state yet; in which case they 

Paper I) JLesaffre et al\ (|2004D, we presented a numeri- ghow & combination of c _ type and j-t VP e features. They 
cal method to compute the time-dependent evolution of migh t also be mildly o r strongly unstable jFim et flil200* 
molecular shocks with a realistic cooling and chemistry m | Srn ith RosenlE o03. and Paper I of this series) , in which 
the presence of a transverse magnetic field. The use of a cage a gteady gtate modd wiU haye limited success 
moving grid algorithm allowed us to reduce the number of 
zones and the computational time. However, the compu- 
tation of the evolution of a stable shock from formation However, different attempts have been made to cir- 
to steady-state still involves one or two days of CPU time cumvent these prob lems. In the fie l d of shocks in su- 
on a 500 MHz workstation. Oscillating shocks require one pernovae remnants, iRavmond et all l|l988l) among oth- 
week or two. This prevents the use of this code to fit shock ers were successful in interpreting spectra with trun- 
parameters to observations. cated steady J-shocks. Those models allowed them to ac- 
On the other hand, steady-state codes give a fast an- count for incomplet e recombination zon es in a filament 
swer (in one or two minutes) and include much richer of the Cygnus loop. IChieze et all l)l998f) discovered that 
physics. They have therefore been extensively used to in- the nascent magnetic precursor in a C-type shock was 

identical to a trunca tion of the steady state sho c k (for 

Send offprint requests to: lesaffre@ast.cam.ac.uk sufficiently late ases) . iFlower & Pineau des Fore"tsl {,1999) 
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used this result to reproduce H2 excitation diagrams in 
Cepheus A West. 

The aim of th i s pap er is to rigorously test the ideas 
oflRavmond et all (|l988h andlFlower fc Pineau des Foretsl 
(1999), by clarifying the relationship between time- 
dependent models and steady-state models. Indeed, we 
will show that time-dependence is within reach of steady 
state codes, as long as the shock is not subject to strong 
instabilities. 

In Sect. 2 we study the stationary equations of mag- 
netic shocks, and derive new analytical relations. Section 
2 is independent of the other sections. Then, in Sect. 3, 
we assess the validity of the stationary approach by using 
the results of our fully hydrodynamical code (Paper I) . In 
Sect. 4, we explain how to build time-dependent models 
of low velocity shocks with the only help of a steady-state 
code. We discuss our results in Sect. 5 and sum up our 
conclusions in Sect. 6. 



2. Analytics of steady shocks 

When the flow is in a steady state, time derivatives can 
be skipped in the steady frame. We derive here a few an- 
alytical relations valid along such steady flows. 



2.1. Dynamical equations in a conservative form 

We recall here the time-dependent monodimensional equa- 
tions of multifluid hydrodynamics with a frozen transverse 
magnetic field. We put them in their conservative form : 



d 



d 



— (rij) + —(rijU n + Jj) = Rj for j neutral specie (1) 
d d 

~Q^( n j) + ~Q~( n j u c + Jj) = Rj f° r 3 ionic specie (2) 
d d 

g^iPnUn) + g^(Pn U n + Pn + n n) = Fc^n (3) 



9 , ^ 9 / 2 B \ r. 



(4) 



8,1 1 



9. , 7 1 



^(^yP„ + ^PnU n ) + Q^M^jPn + ^PnU n + 7T n )] 



= A + Q 1WC + Qi^ c - -u 2 c M e 
1(^ + ^=0 



(7) 
(8) 



n, i, e, and c indices stand for neutrals, ions, electrons, and 
charges, rij, p, u, B, p, and it are respectively the number 
densities, mass densities, velocities, magnetic field, ther- 
mal and viscous pressures. M, F, Q are the mass, momen- 
tum, and heat transfer rates. A denotes radiative losses. 
Rj stands for chemical rates, and Jj for diffusive fluxes. 

2.2. Steady state equations 

Let us assume we are in a frame where the flow is in a 
steady state. We may then drop the Jj terms in equa- 
tions 11181 If we now integrate equations 11181 along the x 
coordinate, we link the state of the gas at one point x 
in the shock to the state of the gas far upstream, i.e. to 
the entrance parameters of the shock (denoted with a 
superscript in the following). We give here the result of 
such an integration in terms of conserved fluxes (dotted 
letters) through the steady region : 
Mass flux : 



M n = p°u° = p u u n - / M n dx' 
Jo 



M c = p° c u° c = Pc u c + / M n dx' 
Jo 

Momentum flux : 



= p 



n ul + p n + 7T n - / F c ^ n dx' 

Jo 

P c = P l { ulf + pl + ^f 



= p c U c + p c + 7Ti 



Energy flux 



8tt 



Jo 



En = <{^L_pl + ip°«) 2 ) 

7 — 1 2 



(9) 
(10) 



(11) 



(12) 



= A n + Qi^ n + Q c ^ n + U n F c ^ n - ~V? n M n (5) 

d, 1 1 2 B 2 . 

^ ( 7-l Pi+ 2^ + 8^ 

d r , 7 1 2 £\i 

+ dx- [Uc{ — Pi+ 2 PiU; + 7Ti + ^ )] 



U n ( TPn + 7iPnU n + 7T n ) - / S n dx' 

7 - 1 2 Jo 



j0\2 



(13) 



TP — vfif 1 r.0 _L ^ rP(nfi\ 2 _L ^ > I 

Fi-u c { -Pi + -Pi(u c ) + — ) 

7 — 1 2 Air 

1 B 2 f x 

= Uc (—L-p i +-p i u 2 c +TTi + - r -)- / Sidx' (14) 
7-1 2 4i j 



1 9 

A; + Q n -^i + Qo^i + u c F n ^ c - -v? c Mi 



d 1 

i-^Pe + §Pe«c) + Q^iM-^jPo + -^PcU 2 c )] 



8 1 
8t y 1 



(6) K = U °(^L- P ° + U(u c ) 2 ) 
7 — 1 z 



= u c (— 1 — p c + l-p c ii 2 + 7T C ) - f S c dx' 
7-12 J 



(15) 
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where S n , S- u and S e stand for the source terms in the 
right hand side of the conservative form of the total energy 
equations 11181 We define A = S n + Si + S c and we note 
that A = A n + A 5 + A e . 
Magnetic flux : 



B = u°B° = u c B 



(16) 



Integrals involve the source terms describing colli- 
sional, chemical, and thermal exchanges between different 
fluids, and radiative losses. The other terms describe the 
conservative phenomena that share mass, momentum, and 
energy between their available reservoirs (thermal, kinetic, 
viscous, magnetic...). 

In each sector of a steady J or C shock, we will now 
get algebraic relations between the dominant conserved 
quantities. We tackle successively the following features, 
in the order in which a parcel of gas entering a magnetised 
shock would meet them : 

— the cold magnetic precursor, in which the ion velocity 
is mainly decelerated, and friction starts to brake the 
neutrals, 

— the hot magnetic precursor, in which the friction has 
brought the temperature to a sufficiently high level 
that H2 cooling starts to play a dominant role (usually 
above neutral temperatures of 10 3 K), 

— the adiabatic front, in which viscosity in the neutrals 
converts their remaining kinetic energy into heat, 

— the relaxation layer, in which the gas cools down, is 
compressed, and gets back to a thermal and chemical 
equilibrium. 

Finally, we derive analytical properties of the atomic 
plateau that follows dissociative shock fronts. 



2.2.1. Cold magnetic precursor 

We assume here and in all the following that the ionisation 
fraction is very low : 



M = M a + M c = p n u n 



(17) 



Viscous, ram and thermal pressure of the charges also 
are negligible compared to neutrals and magnetic pressure. 
Furthermore, since this region is far from the adiabatic 
shock front, neutral viscosity can safely be neglected as 
well. The total momentum flux is then : 



B 

P = P n + P c = p n u 2 n +Pn + ^ 

an 



(18) 



If the temperature of the neutrals stays low, radiative 
losses can be neglected, and the total energy flux is : 

E = E n + E c = Un (^— Pn + \p n ul) + u c f- (19) 

7 — 1 2 4-7T 

Finally, conservation of the magnetic flux through the 
steady region gives : 



B = Bu c 



(20) 



We thus get 4 equations with 5 unknowns p n , u n ,p ni B, 
and u c . One variable can then be chosen to get expressions 
for all the others. For example, p n is solution of a quadratic 
whose coefficients depend on the shock parameters and the 
neutral density p n : 



Pi 



2tt M 2 



7 



B 2 Pi (7 - 1) 



-P 



2ttM 27 
B 2 p n 7 - 1 



M 2 

Pn 



^ ■ 1M 3 

3^- 2d 



= 



(21) 



As another example, u n is solution of a quadratic 
whose coefficients depend on the shock parameters and 



7+1 
2(7-1)' 



Mui 



7-1 



(P 



B 2 

87TU? 



E 



B 2 
4t:u c 



0(22) 



These relations hold up to the point above which ra- 
diative losses cannot be neglected anymore (they would 
hold as well in C-shocks upstream this point). In our sim- 
ulations (Paper I), this corresponds to the point where 
T n > 10 3 K. The gas then enters the hot magnetic precur- 
sor, where H2 cooling becomes dominant. 

2.2.2. Hot magnetic precursor 

In this part of the magnetic precursor, ram pressure is 
directly transferred into radiation via friction. There is 
no more increase in the thermal pressure. Therefore, the 
neutral thermal pressure becomes very quickly negligi- 
ble against magnetic pressure. It also remains negligible 
against ram pressure in the rest of the magnetic precur- 
sor. This leads to the following reduced set of equations : 



M = p n u n 
P 

B = Bu, 



PnK 



Ml 

cStt 



(23) 



We can derive from this a relation between the speeds 
of neutrals and charges that is valid in magnetic precur- 
sors (or C-shocks) downstream the point where | J Adx\ 
dominates over u n —^iP n (usually, when T n is greater than 
10 3 K) : 



M K 8ttu 2 / 



(24) 



This last equation is complementary to eauation l22l These 
equations provide a powerful way to test if the magnetic 
field compression is correctly treated in a multifluid code. 

2.2.3. Adiabatic shock front 

In the shock front, the viscous pressure ir n is one addi- 
tional unknown. We will therefore assume that we know 
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the magnetic field B p at the end of the magnetic precur- 
sor, as well as the amount of energy radiated away in the 
precursor E p — Jp recursor Adx' . Since the shock front is 
very tenuous, it is fair to assume that neither the magnetic 
field nor the integrated radiative losses will vary across it. 

We can then define the new conserved fluxes for this 
region : 



P' 
E' 



P 

E 



Er, 



(25) 



Four equations then combine together : 



M 
P' 

E' 
B 



PnU n 



Pn 



«n(z3iPn + ^Pnul + 7T n ) + U, 



c 8tt 



(26) 



= u c B rj 



Here, we explicitly deduce the pressure in terms of the 
density : 



Pn = (7 



2 Pn 



E' - BB p 
M 



1 \ 

-Pn) 



(27) 



We are not aware of any previous analytic expression 
relating pressure to density throughout an adiabatic shock 
front, even in the absence of magnetic fields. This relation 
is useful to test a code in a shock front. 

In addition, the post-shock velocity u n can be calcu- 
lated by setting 7r n = in equations 1261 which gives the 
following quadratic equation : 



7 P'u n + (7 -Ab p /8tt) = ( 28 ) 



Without magnetic field and energy losses, this quadratic 
gives the post-shock velocity and hence the usual com- 
pression factor in an adiabatic shock. The same kind of 
reasoning will also provide us with analytic predictions in 
dissociative shock fronts (see I2.3|) . 

2.2.4. Relaxation layer 

Here, radiative losses are not negligible anymore, and the 
equation of conservation of energy is left aside. But out- 
side the shock front, viscous pressures are negligible, so 
we only need to make an assumption about the magnetic 
field. Since we neglect the thermal and ram pressure of 
the charges, the momentum conservation of the charges 
yields : 



8tt 



= Pr. 



(29) 



For low shock velocities, the last integral is dominated 
by the magnetic precursor, where most of the ion decelera- 
tion occurs, and B = B p is also a correct approximation in 
the relaxation layer. We arc then left with two equations : 



M = p n U n 

P' = p n ul + p n 



(30) 



Pn(pn) follows easily 



Pn = P' - 



Pn 



(31) 



Note that the intersection of this relation with the ther- 
mal equilibrium relation p n (p n ) gives the final steady post- 
shock conditions. Similarly, the intersection of the alge- 
braic equations in two adjacent sectors of the shock gives 
the physical conditions at the transition between the two 
regions. 

In cases of high shock velocities (greater than 30 kms -1 
for the models in Paper I), a recoupling of the velocities of 
charges and neutrals may happen in the relaxation layer, 
which builds up an additional magnetic pressure. After 
the recoupling zone, we can assume that u n = u c : 



M = p n U a 

P = p n ul + p n 

B = BUn 



Ml 

8tt 



(32) 



Pn{pn) follows easily 



Pn 



p 



Pn 



B 2 



8nM 2 



A 



(33) 



At the high density end of the relaxation layer, the 
total pressure is dominated by magnetic pressure. This 
prevents the gas to be compressed while it cools down, and 
leads to much lower compression factors. The final state of 
the gas in this case corresponds to the steady isothermal 
compression factor calc ulated in the case of mag netic field 
coupled to the gas (see lDraine fc McKeeHl993l equation 
2.19a). 



2.3. Adiabatic dissociative shock 

In an adiabatic dissociative shock, energy losses can still 
be accounted for by conserved quantities since thermal 
energy of the gas is used to dissociate H2, i.e. converted 
into internal energy. We then have : 



E 



u n{^ZZlPn + \pnU n + 7T„) + J** R^Qdx' 



(34) 



where Q = 4.48 eV is the binding energy of the H2 
molecule and R^ 2 is the local dissociation rate. The in- 
tegral term is related to the change in the flux of H2 
molecules by : 



/ Rn 2 Qdx' = Q(0.5 
Jo 



/(*)) 



Mxn 

m 



(35) 



where / is the fractional abundance of H2 by number rel- 
ative to nn (initially / = 0.5), xh is the (fixed) elemental 
hydrogen mass fraction, and pn is the mass of one hydro- 
gen atom. 

For more clarity, we neglect here the effects of the 
charged fluid, but it would be straightforward to include 
the magnetic pressure and radiative effects using primed 
shock parameters like in the previous subsections. We then 
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obtain a modified version of equations 1261 including the 
H2 dissociation energy : 



M 

P 

E 



PnU n 
PnK 



Pn 



Un(-^Pn + kpn< + *„) + (1 - 2/)^ 



(36) 



Unfortunately, we do not have a fixed relation between 
p n and p u as in the non-dissociative case, because the H2 
fraction / varies across the front and adds to the un- 
knowns. But the conditions at the end of the front can 
be found if one sets 7r n = and / = /* where /* is the 
H2 fraction at the end of the shock. So doing, we get an 
equation similar to relation 1251 : 



Ml±l«a- 7 P« n + ( 7 -l)(B-(l-2/*) 



MQx H 
2pn 



) = 0(37) 



We simplify this equation by assuming the high Mach 
number regime, for which P — u°M and E = ^(u°) 2 M : 



7+1 2 , 7-1, on2 n 

— — u n -ju n u +a——(u) =0 



(38) 



2.4. Validation of the analytical results on examples 

We verified the analytical relations derived above by com- 
paring with numerical magnetic shock simulations from 
Paper I. 

First, we find that the temperature and density of the 
atomic plateau are well predicted by the formulae 1401 in 
most dissociative shocks. This was expected, since the 
width of the adiabatic shock is so small that it is very 
likely to be in a steady state : the sound crossing time 
of this feature is much shorter than the time of variation 
of the entrance shock speed. In addition, in strongly os- 
cillating shocks (weakly dissociative case of Paper I) , the 
maximum expansion of the front, corresponding to the 
adiabatic, fully dissociative phase, coincides with the ve- 
locity given by equation 1411 fsee Fig. 4d of Paper I). 

We also verified the relations predicted for non- 
dissociative shocks. As an example, in Fig. ^ we check 
the relations |22] and [2] against the final steady-state of 
a C-type shock (diamonds). The agreement is very good, 
confirming that magnetic compression is correctly treated 
in the code. 



where we defined a = 1 — (1 — 2/*) ( M ?)2 I ^ H ■ 1 — a measures 
the relative decrease in energy flux through the shock front 
due to H2 dissociation. In a non-dissociative shock, a = 1 
since /* = 0.5. a reaches a minimum of 0.74 for a fully 
dissociative shock (/* = 0) just above the dissociation 
limit u° = ltd — 30 kms -1 , and tends to one again for 
very high shock speeds (where the H 2 dissociation energy 
becomes negligible compared to the kinetic flux). 

The compression factor through such a shock can now 
be computed : 



b = 0.1 n=10 4 /cm 3 u=30 km/s 



c 



Pn 



7 +y/7 2 (l-Q)+' 
(7 - l)a 



(39) 



The usual compression C a = is recovered when a = 1 . 
C = 5.8 for 7 = 5/3 and a = 0.74. 

We can get a simple expression for the temperature 
of the atomic plateau (/* = 0) if we neglect the post- 
shock ram pressure (so that p n — p°(u ) 2 ). This is only 
20% accurate for the compression factor obtained (we use 
7 - 5/3) : 



KB 



(40) 



where /i p is the mean molecular weight in the plateau. 
T p is hence nearly quadratic in the entrance velocity for 
strong shocks. 

The knowledge of the compression factor, along with 
the assumption of steadiness of the adiabatic front pro- 
vides us with the velocity of the front relative to the piston 
in the adiabatic phase of a dissociative shock front : 



v = u/(C-l) 



(41) 



0.1 1 

0.1 



Cold magnetic precursor 
Hot magnetic precursor 



1.0 10.0 
Charges velocity (km/s) 



Fig. 1. Steady state analytic relations between the veloc- 
ities of neutrals and charges are compared to an overlaid 
steady C-shock (diamonds). The parameters of the shock 



are u = 20 kms 



10 cm A and b = 0.1, time is 



t = 10° yr. The velocities in the shock frame are com- 
puted using a velocity of the shock front of 0.13 kms -1 , 
inferred from Fig. 3d of Paper I. 



Finally, in figure |21 we plot in diamonds the state of 
the gas for the same shock in a snapshot at t = 100 yrs, 
i.e. well prior to steady (C-type) state. We overlay the 
algebraic relations previously derived for the various shock 
regions. The agreement turns out to be very good. This 
comforts us in the ability of the code to reproduce the 
conservation equations. It also points out the fact that 
steady equations may well be valid to describe a shock at 
early times, even though it has not reached its final steady- 
state. It is to address the domain of validity of this "quasi- 
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steady" approximation that we set up the technique and 
tests described in the next section. 

b=0.1 n = 10Vcm 3 u = 20 km/s 

- \ / 

10" 5" if ? 

E ; o 




Thermochemical equilibrium 
Cold magnetic precursor 
Adiabatic shock from 
Relaxation layer 



10 4 

lO- 2 " 



io-" io-'» 10-" 

Neutral density (g/cm 3 ) 



Fig. 2. Steady state analytic relations between pressure 
and mass density are compared to an overlaid future C- 
shock (diamonds). The parameters of the shock are u = 
20 kms -1 , n = 10 4 cm -3 and b = 0.1, time is t = 10 2 yr. 
The additional necessary parameters B p and E p are read 
in the shock model at the end of the precursor (they are 
not fitted). 



3. Validity of the quasi-steady assumption in 
time-dependent shocks 

In this section we develop a method to characterise the 
local steady velocity for each variable separately, and we 
use it to test the "steadiness" of various shock regions in 
the simulations of Paper I. 

3.1. Local steady velocities and quasi-steady state 

Consider y, one of the TV + 6 state variables that enter the 
set of eauations lllSI Its evolution equation in the frame of 
the piston can be cast in the following form : 

dy_ ^ d(uy) 

at 

or u c is the velocity of the fluid associ- 



(42) 



dx 

where u = u 

ated with y, and s is a term that does not depend on the 
reference frame. 

We define the local steady velocity v y for variable y 
as : the velocity of a reference frame in which the time 
derivative in the evolution equation vanishes. Hence : 



^[{u-v y )y]=s 



When y does not involve a velocity, i.e. does not de- 
pend on the reference frame, a direct expression follows 
for the velocity v y : 



d{uy) 
[ dx 



*]/ 



dy_ 

dx 



(44) 



One has to be more careful when y depends on u. For 
example, in the case y = pu, v y is given implicitly by the 
following quadratic equation : 



y dx 



2V ^ 



dy d{uy) 



dx 



(45) 



The expression 1441 is singular when §f = 0. This can 
easily be understood : if the profile of y is flat, any velocity 
will do. It is therefore crucial to take into account the finite 
numerical precision when trying to evaluate v y with these 
expressions. Indeed, roundoff errors can make the gradient 
of y non zero even in places where it should be. 

Finally, expressions |33| an d EH yield a way of charac- 
terising the "steadiness" of the flow. At a given position 
x, if v y does not depend on y, then there is indeed a frame 
moving at a velocity v(x) in which none of the variables is 
changing in time. Furthermore, if this velocity is constant 
throughout an extended region, then this whole region is 
moving "en-bloc" at velocity v and can be modelled with 
a truncated steady-state model. We say that this region is 
in a quasi-steady state. 

3.2. Validity of the quasi-steady state in 
time-dependent shocks 

For a selection of time steps of each of the dynamical mod- 
els that we simulated in Paper I, we computed the steady 
velocities v y in each zone for each variable y. 

We evaluated the numerical noise in the following way : 
we computed the change Sv yy > in the steady velocity v v 
when each variable y' was changed by 10~ 4 in relative 
value (corresponding to our guess for the numerical preci- 
sion). We then estimated the numerical noise o~ y on vari- 
able y by : 



v' 



{5v yyl ) 



(46) 



The noise- weighted mean steady velocity over all variables 
of a subset S was then computed in each zone, as well as 
the corresponding numerical noise c no i sc : 



(Ev^vcE 1 /^) 



yes 



yes 



yes 



(47) 



(48) 



For charges and neutral momentum, both roots of 
a uadrat ic 1451 where included. But for reasons that will be- 
come clear in the next subsections, magnetic field was ex- 
cluded from this mean. Finally, the variance of individual 
values about this mean velocity and the corresponding 



(43) 

y 2 were estimated 



yes 



( v y 



7^]/(^ 

yes 



V«2) 



1 2/ 



2 

noise 



l G y 



(49) 



(50) 



yes 
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where is the number of variables in the subset S. 

If the numerical noise is well estimated, a value of 
X 2 — 1 indicates that the dispersion of invidual v y val- 
ues is consistent with local numerical noise, i.e. that there 
may exist a common steady velocity for all the variables 
of the subset S at that position. Regions where this is ful- 
filled and v is constant are the quasi-steady regions for the 
set S. 

A less strict criterion can be chosen for the local steadi- 
ness if we think of the ratio vja as a "signal to noise" : 
even if x 2 is high, it may be possible that the ratio vja is 
high. In this case, the steady velocities corresponding to 
different variables are not equal, but they are close to one 
another : therefore, a common steady velocity v is a good 
approximation. 

It turns out that we compute very low values of \ 2 
in quite a few zones. This indicates that the numerical 
precision we have in these zones is far better than our 
estimate of 10 -4 . We hence rather use the criterion based 
on the ratio v j a and define the quasi-steady regions as the 
regions with a constant velocity v and a good "signal to 
noise" ratio. 

The subset of variables S on which the averages are 
computed should be the whole set of variables. However, 
the last section of this paper needs only that the dy- 
namically important variables be in a quasi-steady state. 
Therefore we present here the results for a subset S includ- 
ing the temperatures of the three fluids, the four velocities 
(roots of equation 1451 for y = p n u n and y — p c u c ), He, 
H2, H2O, CO and OH densities. We do not include the 
magnetic field in magnetised shocks, because its steady 
velocity differs from the other variables as we will show 
in the next section. We also did the calculation for S in- 
cluding all the variables but the magnetic field, and found 
that it did not change the general conclusions : the signal 
to noise ratio is slightly less good, and a few zones are 
not quasi-steady anymore because of marginal chemical 
species having a different steady velocity than the bulk of 
the variables. 

We summarise the results of our investigation in 
the next two subsections, devoted respectively to non- 
magnetised and magnetised shocks. 

3.2.1. J-type shocks (-5=0) 

N on- dissociative J-type shocks 

Figure |3| shows the x 2 as well as v with error bars ±<r 
in each zone of a typical snapshot of a non-dissociative J- 
type shock. The adiabatic front and relaxation layer show 
a good "signal to noise" ratio and a flat steady velocity. On 
the contrary, the pre-shock zones show a huge dispersion 
around the steady velocity, expected from the homogene- 
ity of the medium there. The fact that we retained both 
velocities from expressionl45ldoes not alter the results, be- 
cause the numerical errors on velocities are actually much 
larger than for the other variables. 











Adiabatic front Relaxation layer 



20 40 60 ao 100 

Zone number 

Fig. 3. We plot the x 2 values and steady velocities for 
each variables of S in each zone of snapshot t — 200 years 
for the shock with parameters b = 0, n — 10 4 cm -3 , and 
u = 20 kms -1 . Each dot represents the steady velocity of 
one variable computed thanks to expressions ^] and 1451 
The error bars are via evaluated zone per zone on these 
values. We indicate the computational domains associated 
to the adiabatic front and the relaxation layer of the shock. 

We find that all of our non-dissociative J-type shocks 
are in a quasi-steady state from the adiabatic phase to 
the steady phase (only the initial formation of the shock 
front is not quasi-steady). Therefore, a snapshot of such a 
shock will always coincide with the truncated structure of 
a J-type steady state. 

In Paper I, we plotted the trajectory of the point of 
maximum ratio of viscous over thermal pressure (see Fig. 
lb). We compute here the average steady velocity v over 
the whole structure of the shock (trimmed from the pre- 
shock values) at various times and overlay it over this 
trajectory in Fig.0] Error bars show the good consistency 
of the test, and the correspondence of v with the velocity 
at which the viscous shock front moves away from the 
piston. 

Note that the entrance velocity in the shock is not the 
upstream velocity u of the fluid towards the piston, but 
rather vP = u + v. Hence, the entrance shock speed for the 
truncated steady shock is evolving in time. In Sect. 01 we 
present a way to reconstruct this evolution. 

Dissociative J-type shocks 

On the contrary, dissociative shocks are almost never 
in a quasi-steady state. For weakly dissociative velocities, 
we could not come up with a coherent picture. This was 
expected, since these shocks are highly unstable with large 
bouncing oscillations between a fully dissociative expan- 
sion phase and a non-dissociative recoil phase (Paper I). 
Partly ionised shocks are less unstable. We illustrate their 
behaviour with a typical example shown in Fig. [5j The 
adiabatic front is generally in a quasi-steady state with 
velocity roughly equal to the velocity of the viscous max- 
imum, but the dispersion is much higher than for non- 
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Fig. 4. Trajectory and velocity away from the piston of 
the J-shock with parameters b — 0, n — 10 4 cm~ 3 , and 
u = 20 kms -1 (from Paper I). Overlaid diamonds are 
the steady velocities v averaged over all variables and all 
zones, for each snapshot analysed. 



dissociative shocks. The first plateau that follows (with 
H ionisation and Lyman cooling) is not at all in a quasi- 
steady state : the dispersion is huge and the mean velocity 
is not even constant. The second plateau (H recombina- 
tion) seems rather quasi-steady, but with a velocity much 
lower than the adiabatic shock front. This warns us that 
steady-state diagnostics may be hopeless for weakly non- 
dissociative shocks, and that we have to be cautious for 
partly ionising shocks. 



3.3. CJ-type and C-type shocks 

Figure [B] shows v with error bars ±cr in each zone of 
an early snapshot of a non-dissociative C-type shock. 
The shock then has a composite CJ structure made of 
a magnetic precursor, followed by a non-dissociative adi- 
abatic front and a relaxation layer. We find that the mag- 
netic precursor is in a quasi-steady state with a high 
steady velocity. Th is is in agreement with the remark of 
IChieze et all ( 1998t) . At the very end of the magnetic pre- 
cursor, the large dispersion in steady velocities is due to 
the fact that two different steady velocities coexist among 
the variables. Following this, the adiabatic front and re- 
laxation layer appear to be in a quasi-steady state, but 
with a much lower velocity than the magnetic precursor. 
However, strictly speaking, it is not a real quasi-steady 
state, as the steady velocity for the magnetic field remains 
equal to the high steady velocity of the magnetic precur- 
sor. 
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Fig. 5. Same as Fig. [31 but for a partly ionising shock of 
parameters b — 0, n — 10 4 cm~ 3 , and u = 40 kms -1 at 
time t — 220 years. Here, the scale of the plot is linear, 
so that the dispersion is in fact much greater than for the 
non-dissociative shocks, even in the second plateau. 



Fig. 6. Same as Fig. |3| but for time t = 100 years of 
the future C-shock with b = 0.1, n = 10 4 cm 



-3 



u = 20 kms 1 . 
magnetic field. 



and 

We also show the steady velocity for the 



Figure [7| shows that the two steady velocities (relax- 
ation layer and precursor) correspond well to the velocities 
of the viscous maxima of neutrals and charges determined 
by time-derivation of their trajectory. As time evolves, the 
velocity of the magnetic precursor and that of the relax- 
ation layer get closer to one another, and finally coincide 
after the J-front has disappeared. The C-type structure is 
then in a quasi-steady state as a whole. 

In principle, one should then be able to model an early 
age of a low velocity C shock by combining a truncated 
C-type model with a truncated J-type model (in which 
the magnetic field is treated appropriately). The problem 
is a bit more complex than in the J-type case because we 
need now to determine two different truncation distances 
and two sets of entrance parameters, but we will show in 
Sect. l4.2l that it is possible to solve. 
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This picture holds for all shocks with magnetic field, 
as long as there is no dissociation or ionisation plateau. 
In the case of dissociative velocities, the same problems 
described in the previous subsection arise in the corre- 
sponding features of the relaxation layer. 




100.00 

Time (years) 



Fig 

b = 



7. Same as Fig. 0J but for a C-shock of parameters 



0.1, n 



10 cm and u 



20 kms . Curves 



plot the trajectory and velocity of the neutral (solid) 
and charged (dashed) viscous fronts. Diamonds are the 
quasi-steady velocities of the relaxation layer and adia- 
batic front. Triangles are the quasi-steady velocities of the 
magnetic precursor. 



one would then like a method to relate these parameters 
to the parameters of the shock in the piston frame, and to 
the age of the shock. Conversely, one would like to build 
at will a snapshot of a shock of given age and parameters 
(in the piston frame), with the only help of a steady state 
code. We come up with such a procedure in the following. 

A steady-state code provides us with the steady profile 
of any variable in the frame of the shock front for a given 
set of entrance parameters (u , n). Say, the steady velocity 
u s (u°,n;x), where x is the distance from the shock front. 
If we are given an inflow speed and density (it, n) in the 
frame of the piston and a time t, the problem is to find 
what is the entrance velocity u° at the same time t in 
the frame of the shock front as well as the corresponding 
distance r between the shock front and the piston. 

To help set up the notations in both frames of the 
shock and of the piston, we sketched in Fig.|Hlthe different 
lengths and velocities involved. 



Piston frame 



Shock frame 



UpJ = 



r Or 

Fig. 8. Schematical view of a J-type shock in the piston 
frame and in the shock frame. 



This picture is the same for all magnetic shocks. The 
only difference between CJ-type and C-type is whether or 
not the J-front has disappeared when the steady veloc- 
ities of the relaxation layer and the magnetic precursor 
converge. From this remark, we will obtain in Sect. POTl 
a way to assess if a low velocity shock will eventually be- 
come a steady CJ-type shock. 

4. Time-dependent constructions of shocks at 
early times 

Here, we derive methods of reconstruction of time- 
dependent shocks using truncated steady models. Those 
constructions will be meaningful only for the shocks in 
which the quasi-steady state has been validated at all 
times, although they can in principle be realised in any 
shock. Due to their different complexities, we treat suc- 
cessively the case of non-dissociative J-type shocks and 
non-dissociative magnetised shocks. 

4.1. Non-dissociative J-type shocks 

Section \'6.2. II has shown that the whole structure of non- 
dissociative J-type shocks is at all times quasi-steady. One 
may then safely fit truncated steady models to observa- 
tions. The fitted parameters would be the entrance param- 
eters in the shock frame and the truncation distance. But 



At any given time, velocities in the shock frame are 
found by adding v = r to velocities in the piston frame. 
The entrance shock speed is then : 

u° = u + r, (51) 

while the velocity at the piston — which must be null in 
the piston frame — is, in the J-shock frame : 

u s (u°, n; r) = r. (52) 

These relations combine to give an implicit equation link- 
ing r to r : 

r = u s (u + r,n;r) (53) 

Furthermore, in a quasi-steady state, mass conserva- 
tion requires that u s (u°,n;r) x C = u , i.e. r = v = 
u/(C — 1) where C is the compression factor at the pis- 
ton. From the adiabatic phase (C = 4) to the steady state 
(C « 1), the speed v of the front thus decreases from 
■|u to nearly 0. Therefore, if the steady state code pro- 
vides u s (u° 7 n;r) for a range of velocities u < u° < |it, 
the equation 1531 is an implicit ordinary differential equa- 
tion straightforward to integrate up to time t with initial 
conditions r = and r = Actually, an interpolation 
between a few steady models might be sufficient to get 
accurate results. 
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Conversely, if the problem is to recover the time from 
a steady model truncated at distance R, we only have to 
integrate equation 1531 backward in time up to the point 
where r = and compute t — A dr. 

High compression factor approximation : 
For high compression factors, the final f is small 
enough to be neglected with respect to u, which makes the 
computation even easier. The integral yielding the age is 
dominated by the very low velocities, which are also the 
most recent ones. Therefore, we only need one steady state 
model u s (x) = u s (u, n, B; x). The age of such a truncated 
shock is simplified in the following way. 

t= / — dx (54) 
Jo u s 

In this last expression, one recognises the flow time across 
the shock. Since for strongly radiative shocks, the com- 
pression factor rises very quickly, this approximation is 
valid even for very young ages. For example, the shock of 
parameters b = 0, n = 10 4 cm' 3 , and u = 20 kms -1 has 
a compression factor of 10 at as early as t = 1 year (see 

Fig. mi . 

4.2. Non-dissociative magnetised shocks 



Piston frame 




X 



Fig. 9. Schematical view of an early magnetised shock in 
the piston frame. 

The analysis of Sect. EH showed that low- velocity mag- 
netised shocks are composed of two quasi-steady regions : 
a magnetic precursor, and a non dissociative J-type fea- 
ture. In faster shocks, where the entrance velocity in the 
J-type feature is dissociative, the J-type structure is not 
in a quasi-steady state, although the magnetic precursor 
is. We thus restrict our analysis to the non-dissociative 
cases. 

In principle, one should then be able to model a snap- 
shot of such a shock by gluing together two truncated 
steady C and J models. The problem is to determine the 



entrance parameters and lengths of each of the two shock 
features, for a given time t and a given set of parame- 
ters (u, n, B) in the piston frame. A rigorous construction 
method is outlined below. 

In the following, rc and rj denote respectively the dis- 
tance of the fronts of the C-type and the J-type features 
with respect to the piston (see Fig.|5J|. Variables computed 
in the reference frame of the C-type feature are denoted 
with C subscripts. Variables in the J-type feature frame 
are specified with J subscripts. Entrance parameters in 
the shocks have a superscript. 

Now, let us assume that we know the positions rc and 
rj at current time t. To solve for their evolution, we need 
to find equations that will determine rc and fj. 

The entrance velocities in the C-type feature are sim- 
ply: 

u "c = r c +u 

n c — n 
B° c = B 

The steady-state code for the C-type feature then provides 
us with the entrance values of velocities, densities, and 
magnetic field at the position of the J-front, x — rc — rj. 
After a suitable change of reference frame for the veloci- 
ties, they determine the entrance parameters in the J-type 
shock : 

w°j = rj - r c + u nC (u + r c ,n, B; r c - rj) 

"cj = O - r c + u cC {u + fc,n, B; r c - rj) 

n °J = n nC (u + r c ,n, B;r c - rj) (56) 

n°j = n cC (u + fc,n, B; r c - rj) 

Bj = B c (u + rc, n, B; r c - rj) 

The steady-state J-shock must then be integrated. A 
multifluid treatment is necessary, since u n and u c at the 
entrance of the J-front are different. Furthermore, a spe- 
cial treatment of magnetic field compression is necessary, 
since Sect. 13.31 showed that the steady velocity for the 
magnetic field in the J-type feature is not rj, but remains 
the same as in the magnetic precursor, namely vb = rc- It 
means that the product of B with the velocity of charges 
computed in the frame of the C shock remains constant 
through the J-type feature : 

Bx(u c j-fj + r c ) = B c xu cC (57) 

Therefore, the evolution of the J-type feature depends not 
only on the entrance parameters determined above but 
also on fc — r j . 

As in the non-magnetic case, the derivatives rc and 
fj are then determined by stating that both charges and 
neutrals have to be at rest near the piston, i.e., in the 
frame of the J-front : 

rj = UnjKj.uSjj^nSj.BS.rc-rjjrj) ( , 

I d0 • • \°°) 

We thus get two independent implicit equations for fj 
and rc . Numerical techniques to solve these equations still 
need to be designed, but should not be too hungry in CPU 
time. 
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4.2.1. Final steady state: C or CJ ? 

We have now a means of computing the time evolution 
of a magnetised shock with only a steady-state code. One 
should then be able to integrate it from initial conditions 
rj = tq = up to the steady state where fj and fc are 
equal constants. If during the evolution the entrance veloc- 
ity in the J-shock u°j becomes subsonic, then one should 
stop the integration because the J-shock disappears, and 
the remaining sound wave propagates through the struc- 
ture until a stationary C-type structure is obtained. If j 
stays supersonic when f j and re are equal constants, then 
the result is a CJ-type shock steady-state. Hence, one is 
forced to integrate over time r'j and r'c given implicitly 
by equations 1581 to know what is the final steady-state 
corresponding to a given set of parameters (it, n, B). 

However, one could also think about solving equations 
1561 and 1581 for given arbitrary values of the final front ve- 
locity Vf = f j = fc . The result would be a series of physi- 
cally consistent steady CJ-type states, each characterised 
by a different distance re — rj. Only one of these is se- 
lected by the time evolution but, if the entrance param- 
eters (u, n, B) are allowed to evolve in time, it might be 
possible that several (or even all) of these final states can 
be realised. The final state would then depend on the evo- 
lution history of the entrance parameters. 

In fact, a very easy way to exhibit one of those CJ-type 
steady-states would be to use a multifluid steady-state 
code, and trigger the viscous dissipation in the neutral at 
a given position rc — rj where the neutral velocity is still 
supersonic. 

4.2.2. Low velocity, high compression factor 
approximations 

For all the low velocity cases encountered in our simula- 
tions, we noted that after one year of time, the velocity of 
the charges was already almost brought to rest at the end 
of the magnetic precursor. This approximation yields the 
following equation : 

fc = u c c(u + fc,n,B;rc - rj) (59) 

which implicitly gives the velocity fc- The velocity fj can 
then be retrieved by solving the first equation of the set 

EE 

Just like for the J-type shocks, high magnetic com- 
pression factors will lead to fc negligible before u, and 
will facilitate the integration of equation 1591 In this case, 
and if in addition rj << re, the age of the shock is given 
by: 

/•rc-rj i 

t= / da; (60) 

Jo UcC 

which is the flow time of the charges across the magnetic 
precursor. 



5. Discussion 

The analytic relations we found make good benchmarks 
for testing codes. In addition, they might provide some 
theoretical basis for further investigation of the properties 
of these shocks in the parameter space. 

The quasi-steady state analysis of shocks opens a new 
field of possibilities for the steady-state codes. We compare 
hereafter our method to previously used algorithms, and 
sketch possible extensions of our method. 

5.1. Comparison with previous work 

Our quasi-steady state analysis of J-sho cks justifies the 
use of truncated steady-state J-shocks bv lRavmond et al\ 
l|l988|) . We provide more theoretical basis to link the true 
age of the shock to the truncation distance used. 

iFlower fc Pineau des Foretd l|l999|) and 

iLe Bourlot et all l)2002j) use simple algorithms to 
produce mixed C-type and J-type fe atures to mimic 
time-d ependent magnetised shocks. ILe Bourlot et all 
(120021) greatly improved the method used by 
IFlower &: Pineau des Foretsl l|l999|) since they keep 
the multifluid treatment of the flow through the relax- 
ation layer. They just switch on viscosity in the neutral 
fluid when they encounter a sonic point. The present 
analysis gives a less heuristic way to know at which point 
the viscosity should be switched on, and Paper I has 
already shown that it can be way upstream a sonic point. 
Furthermore, we specify that a change of velocity frame 
has to be done at the end of the magnetic precursor, 
except for the magnetic field equation. Finally, we state 
where the J-type structure has to be truncated for a 
given time t. 

Our new method should therefore lead to more accu- 
rate results, and will allow the construction of much ear- 
lier phases of magnet ised shocks. It sh o ws as well that 
the criterion used by ILe Bourlot et all (|2002l) to assess 
whether steady-states will be of CJ-type (occurrence of 
a sonic point in the neutral fluid) has to be revised. CJ- 
type steady states may in fact occur at lower speeds, when 
velocity recoupling between neutral and charges enhances 
magnetic compression near the piston, and slows down the 
precursor to the expansion speed of the J-front. 

5.2. Possible extensions of the method 

First, let us point out that the time-dependent construc- 
tion method derived here relies only on the quasi-steady 
state assumption for a limited number of variables, namely 
velocities, densities, and magnetic field. For example, if a 
set of chemical species can be identified to have no impact 
on the dynamics, they can be skipped in the process of 
building the truncation radii, and computed only in the 
last resort. 

Following the same idea, if non-dynamically important 
species happen to be non quasi-steady, they can be post- 
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processed in parallel to the quasi-steady time-evolution 
with a Lagrangian code. 

Here we present an algorithm for two kinematic flows 
(charges and neutrals), but the same method may be 
implemented for more flows. Especially, the treatment 
of charged grains could be envisaged, in relatio n to 
the questions raised by ICiqlek fe Robergel l)2002f) and 
iFlower fe Pineau des Foretsl l|2003j) . The only caveat is 
that we do not yet have a consistent check for the validity 
of the quasi-steady assumption. 

Finally, our algorithm is straightforward to apply with 
slowly changing input conditions (u, n, B) in the piston 
frame. One has only to bear in mind that if these pa- 
rameters change over time-scales much shorter than the 
crossing time scale of the shock, then the quasi-steady as- 
sumption is very likely to be violated. 

5.3. Limitations of our method 

Our algorithm is based upon the quasi-steady state as- 
sumption. However, it will give results with any shock. 
One problem is that we still have no other way to assess 
the validity of the steady-state assumption than comput- 
ing the time-dependent evolution with a fully hydrody- 
namical code. 

We encountered several cases where this assump- 
tion was not realised. Strongly unstable shocks like the 
weakly dissociative ones violate strongly this assump- 
tion. Fortunately, they seem to happen for a very re- 
stricted range of parameters. Partly ionising shocks are 
very slightly unstable, and are closer to meet the quasi- 
steady state assumption. They might therefore be ac- 
counted for by our algorithm. Finally, quite a few mag- 
netised shocks have unstable entrance velocities for the J- 
shock only at early times, and are afterwards quasi-steady 
at all times. These shocks may be as well within reach of 
our algorithm if one is ready to skip the early evolution. 
However, one should always be cautious when a plateau 
with dissociated molecules appears in a steady-state com- 
putation. 

An other situation where the quasi-stationary assump- 
tion may be strongly violated is the case where a dynam- 
ically important chemical specie is not in a quasi-steady 
state. This might happen when a dominating cooling agent 
varies on very short time scales. Furthermore, diffusion ef- 
fects, if they turn out to be important, will destroy the 
quasi-steady state as well. 

6. Conclusions 

In a companion paper (Paper I), we produced fully time- 
dependent numerical simulations of molecular shocks. 

In the present paper, we derived new analytical rela- 
tions valid at quasi-steady state, and successfully checked 
them on our simulations. These relations provide useful 
benchmarks to test existing and future multifluid codes. 

In light of the simulations run in Paper I, we investi- 
gated carefully the validity of the quasi-steady state ap- 



proximation. It was found that at all times stable shocks 
could be accounted for by truncated steady models. We 
point out as well that caution has to be kept regarding 
the use of steady-state models for dissociative velocities. 

Finally, we produced a new algorithm based on the 
quasi-steady state assumption. With only a steady-state 
code, this method is able to compute time-dependent 
snapshots of shocks in the presence or not of a mag- 
netic field. Therefore, it brings time-dependence within 
the reach of steady models, and should greatly improve the 
diagnostics of observed molecular shocks. Furthermore, it 
provides a way of assessing the CJ nature of a magnetised 
shock. Finally, this algorithm can be extended to many 
shocks other than molecular, provided that the quasi- 
steady state approximation is validated. 
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